
options(digits=4, width=70)

# make sure packages are installed prior to loading them
library(PerformanceAnalytics)
library(zoo)
library(boot)
library(tseries)

# get monthly adjusted closing price data on SP500, NIKKEI225 and DAX from Yahoo
# using the tseries function get.hist.quote(). Set sample to Sept 2005 through
# Sep 2010. Note: if you are not careful with the start and end dates
# or if you set the retclass to "ts" then results might look weird

startDate="2012-08-12"
endDate="2014-03-12"

# get the last five years of monthly adjusted closing prices from Yahoo!
SP500.prices = get.hist.quote(instrument="%5EGSPC", start=startDate,
                              end=endDate, quote="AdjClose",
                              provider="yahoo", origin="1970-01-01",
                              compression="m", retclass="zoo")
# change class of time index to yearmon which is appropriate for monthly data
# index() and as.yearmon() are functions in the zoo package 
#                             
index(SP500.prices) = as.yearmon(index(SP500.prices))

class(SP500.prices)
colnames(SP500.prices)
start(SP500.prices)
end(SP500.prices)

NIKKEI225.prices = get.hist.quote(instrument="%5EN225", start=startDate,
                              end=endDate, quote="AdjClose",
                              provider="yahoo", origin="1970-01-01",
                              compression="m", retclass="zoo")
index(NIKKEI225.prices) = as.yearmon(index(NIKKEI225.prices))

DAX.prices = get.hist.quote(instrument="^GDAXI", start=startDate,
                             end=endDate, quote="AdjClose",
                             provider="yahoo", origin="1970-01-01",
                             compression="m", retclass="zoo")
index(DAX.prices) = as.yearmon(index(DAX.prices))

RTS.prices = get.hist.quote(instrument="RTS.RS", start=startDate,
                            end=endDate, quote="AdjClose",
                            provider="yahoo", origin="1970-01-01",
                            compression="m", retclass="zoo")
index(RTS.prices) = as.yearmon(index(RTS.prices))


CH.prices = get.hist.quote(instrument="000001.SS", start=startDate,
                            end=endDate, quote="AdjClose",
                            provider="yahoo", origin="1970-01-01",
                            compression="m", retclass="zoo")
index(CH.prices) = as.yearmon(index(CH.prices))


HK.prices = get.hist.quote(instrument="%5EHSI", start=startDate,
                           end=endDate, quote="AdjClose",
                           provider="yahoo", origin="1970-01-01",
                           compression="m", retclass="zoo")
index(HK.prices) = as.yearmon(index(HK.prices))


# create merged price data
lab4Prices.z = merge(SP500.prices, NIKKEI225.prices, DAX.prices, 
                     RTS.prices, CH.prices, HK.prices)

# rename columns
colnames(lab4Prices.z) = c("SP500", "NIKKEI225", "DAX", "RTS", "CH", "HK")

# calculate cc returns as difference in log prices
lab4Returns.z = diff(log(lab4Prices.z))

#
# 3. Create timePlots of data
#
plot(lab4Returns.z, plot.type="single", lty=1:3, col=1:6, lwd=2)
legend(x="bottomleft", legend=colnames(lab4Returns.z), lty=1:3, col=1:6, lwd=2)
abline(h=0)
title("Monthly cc returns")

#
# 4. Create matrix of return data and compute pairwise scatterplots
#

ret.mat = coredata(lab4Returns.z)
colnames(ret.mat)
head(ret.mat)
SP500 = ret.mat[,"SP500"]
NIKKEI225 = ret.mat[,"NIKKEI225"]
DAX = ret.mat[,"DAX"]
RTS = ret.mat[,"RTS"]
CH = ret.mat[,"CH"]
HK = ret.mat[,"HK"]
pairs(ret.mat, col="blue")

#
# 5. Compute estimates of CER model parameters
#
muhat.vals = apply(ret.mat, 2, mean)
muhat.vals
sigma2hat.vals = apply(ret.mat, 2, var)
sigma2hat.vals
sigmahat.vals = apply(ret.mat, 2, sd)
sigmahat.vals
cov.mat = var(ret.mat)
cov.mat
cor.mat = cor(ret.mat)
cor.mat
covhat.vals = cov.mat[lower.tri(cov.mat)]
rhohat.vals = cor.mat[lower.tri(cor.mat)]
names(covhat.vals) <- names(rhohat.vals) <- 
  c("SP500,NIKKEI225","SP500,DAX","SP500,RTS", "SP500,CH",
    "NIKKEI225,DAX", "NIKKEI225,RTS", "NIKKEI225,CH", "DAX,RTS", "DAX,CH", "RTS,CH")
covhat.vals
rhohat.vals

# summarize the CER model estimates
cbind(muhat.vals,sigma2hat.vals,sigmahat.vals)
cbind(covhat.vals,rhohat.vals)

# plot mean vs. sd values
plot(sigmahat.vals, muhat.vals, pch=1:5, cex=2, col=1:3, 
     ylab = "mean", xlab="sd (risk)")
abline(h=0)     
legend(x="topright", legend=names(muhat.vals), pch=1:5, col=1:3, cex=1.5) 


#....
#
# 11. 24 month rolling estimates of mu and sd
#

# rolling analysis for SP500
roll.mu.SP500 = rollapply(lab4Returns.z[,"SP500"],
                          FUN=mean, width = 24, align="right")

roll.sd.SP500 = rollapply(lab4Returns.z[,"SP500"],
                          FUN=sd, width = 24, align="right")

plot(merge(roll.mu.SP500,roll.sd.SP500,lab4Returns.z[,"SP500"]), plot.type="single",
     main="24-month rolling means and sds for SP500", ylab="Percent per month",
     col=c("blue","red","black"), lwd=2)	
abline(h=0)
legend(x="bottomleft", legend=c("Rolling means","Rolling sds", "SP500 returns"), 
       col=c("blue","red","black"), lwd=2)

# rolling analysis for NIKKEI225
roll.mu.NIKKEI225 = rollapply(lab4Returns.z[,"NIKKEI225"],
                          FUN=mean, width = 24, align="right")

roll.sd.NIKKEI225 = rollapply(lab4Returns.z[,"NIKKEI225"],
                          FUN=sd, width = 24, align="right")

plot(merge(roll.mu.NIKKEI225,roll.sd.NIKKEI225,lab4Returns.z[,"NIKKEI225"]), plot.type="single",
     main="24-month rolling means and sds for NIKKEI225", ylab="Percent per month",
     col=c("blue","red","black"), lwd=2)	
abline(h=0)
legend(x="bottomleft", legend=c("Rolling means","Rolling sds", "NIKKEI225 returns"), 
       col=c("blue","red","black"), lwd=2)

# rolling analysis for DAX
roll.mu.DAX = rollapply(lab4Returns.z[,"DAX"],
                          FUN=mean, width = 24, align="right")

roll.sd.DAX = rollapply(lab4Returns.z[,"DAX"],
                          FUN=sd, width = 24, align="right")

plot(merge(roll.mu.DAX,roll.sd.DAX,lab4Returns.z[,"DAX"]), plot.type="single",
     main="24-month rolling means and sds for DAX", ylab="Percent per month",
     col=c("blue","red","black"), lwd=2)	
abline(h=0)
legend(x="bottomleft", legend=c("Rolling means","Rolling sds", "DAX returns"), 
       col=c("blue","red","black"), lwd=2)

# rolling analysis for RTS
roll.mu.RTS = rollapply(lab4Returns.z[,"RTS"],
                          FUN=mean, width = 24, align="right")

roll.sd.RTS = rollapply(lab4Returns.z[,"RTS"],
                          FUN=sd, width = 24, align="right")

plot(merge(roll.mu.RTS,roll.sd.RTS,lab4Returns.z[,"RTS"]), plot.type="single",
     main="24-month rolling means and sds for RTS", ylab="Percent per month",
     col=c("blue","red","black"), lwd=2)	
abline(h=0)
legend(x="bottomleft", legend=c("Rolling means","Rolling sds", "RTS returns"), 
       col=c("blue","red","black"), lwd=2)

# rolling analysis for CH
roll.mu.CH = rollapply(lab4Returns.z[,"CH"],
                        FUN=mean, width = 24, align="right")

roll.sd.CH = rollapply(lab4Returns.z[,"CH"],
                        FUN=sd, width = 24, align="right")

plot(merge(roll.mu.CH,roll.sd.CH,lab4Returns.z[,"CH"]), plot.type="single",
     main="24-month rolling means and sds for CH", ylab="Percent per month",
     col=c("blue","red","black"), lwd=2)  
abline(h=0)
legend(x="bottomleft", legend=c("Rolling means","Rolling sds", "CH returns"), 
       col=c("blue","red","black"), lwd=2)